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Thin membranes exhibit complex responses to external forces or geometrical constraints. A 
familiar example is the wrinkling, exhibited by human skin, plant leaves, and fabrics, resulting from 
the relative ease of bending versus stretching. Here, we study the wrinkling of graphene, the thinnest 
and stiff est known membrane, deposited on a silica substrate decorated with silica nanoparticles. 
At small nanoparticle density monolayer graphene adheres to the substrate, detached only in small 
regions around the nanoparticles. With increasing nanoparticle density, we observe the formation 
of wrinkles which connect nanoparticles. Above a critical nanoparticle density, the wrinkles form 
a percolating network through the sample. As the graphene membrane is made thicker, global 
delamination from the substrate is observed. The observations can be well understood within 
a continuum elastic model and have important implications for strain-engineering the electronic 
properties of graphene. 



I. INTRODUCTION 

Thin films supported on substrates are of technological 
importance and are commonplace in biological systems 
such as cell walls and hard skins on soft plant and animal 
tissues. As the thickness t of a slab of material is reduced, 
it becomes more susceptible to out-of-plane deformation 
(bending) compared to in-plane deformation (stretching), 
resulting in morphological transitions. For example, thin 
films deposited on soft compliant substrates display wrin- 
kling patterns under compressive stress PQ, and mem- 
branes resting on fluids wrinkle by capillary forces [2]. 
This wrinkling is a ubiquitous phenomenon found in sys- 
tems ranging from human skin to draping fabric [3H5] and 
has also been exploited to fabricate flexible electronic de- 
vices urn?]- 

Due to the inextensible but bendable nature of the 
sp 2 carbon bond, graphene's effective mechanical thick- 
ness t e ff = ( YIkJE^d ) 1 ^ 2 is less than 1 A [8], where k « 
1 eV |S] is the bending rigidity and E 2 d ~ 2.12 x 10 3 
eV/nm 2 [10] is the tensile rigidity; this is likely the small- 
est mechanical thickness ever achieved for any material. 
Graphene has been anticipated to exhibit a rich vari- 
ety of wrinkling and delamination behaviors [TTHT6] . Al- 
though graphene adheres conformally to smooth nano- 
scale features with high fidelity [17], graphene wrinkling 
has been observed under compressive [18 or tensile [19] 
stress caused by thermal cycling, and graphene on period- 
ically corrugated elastic [20] and metal [21 substrates ex- 
hibit transitions from adhesion to delamination. Morpho- 
logical features such as wrinkles and conical singularities 
may produce non-uniform strain in graphene [22] , which 
produce both scalar and vector potentials and mimic the 
effect of a magnetic field on graphene's electronic struc- 



ture [23] . Recent experimental results suggest these effec- 
tive "pseudomagnetic" fields can exceed several hundred 
Tesla [24 , and graphene devices based on "strain engi- 
neering" have been proposed [25 , 26]. Understanding the 
mechanical response of graphene to non-uniform stress is 
thus a critical first step toward strain engineering its elec- 
tronic properties. 

Here we report a systematic study of the wrinkling 
of graphene membranes supported on Si02 substrates 
with randomly placed topographic perturbations, pro- 
duced by Si02 nanoparticles. The wrinkling is probed 
as a function of nanoparticle density p np and membrane 
thickness (using multilayer graphene). At low p npi mono- 
layer graphene largely conforms to the substrate ex- 
cept for small regions around the nanoparticles, where 
graphene is detached. Wrinkles (or folds or ridges) form 
as p np increases, connecting pairs of protrusions. The 
observed maximum wrinkle length is predicted quantita- 
tively within a simple elastic model. Above a critical den- 
sity the wrinkles percolate to form a network spanning 
the entire sample. As the thickness of graphene increases, 
it stiffens and delaminates instead of wrinkling. Since the 
wrinkling acts to remove inhomogeneous in-plane elas- 
tic strains through out-of-plane buckling, our results can 
be used to place limits on the possible in-plane strain 
magnitudes that may be created in graphene to realize 
strain-engineered electronic structures [22H26] , 

II. EXPERIMENTAL DETAILS 

Silica-nanoparticle colloidal dispersions (Nissan Chem- 
ical America Corp., Snowtex-O) were diluted to various 
concentrations of 0.5—3.0 wt % by deionized water (Fisher 



FIG. 1. AFM images (lxl /xm 2 ) of graphene on SiCb nanoparticle/SiCb substrates for a nanoparticle density of (a) 11, (b) 
22, (c) 49, (d) 90, and (e) 170 /im" 2 , respectively. 



Scientific, Water HPLC Grade). The diluted suspensions 
were sonicated for 30 min in a water bath to break ag- 
glomerations before spin-coating the nanoparticles onto 
a substrate. Spin-coating was performed on Si substrates 
with a 300 nm-thick oxide layer at 4000 rpm for 30 sec- 
onds. The density of nanoparticles on substrates ranged 
from 2 /im -2 to 258 /im -2 , depending on the concentra- 
tions of the nanoparticle dispersions. After spin-coating, 
the samples were completely dried on a hotplate at ~ 
150 °C for 2 hours. Graphene flakes were mechanically 
exfoliated from Kish graphite onto Si02 substrates cov- 
ered with the silica nanoparticles (mean diameter 7.4 ± 
2.2 nm; Ref. 27). Thicknesses of graphene films were 
identified with an optical microscope, atomic force mi- 
croscopy (AFM), and/or Raman spectroscopy. The sizes 
of graphene sheets were typically more than 10 /im x 
10 /xm, which were much larger than an estimated dis- 
tance between nanoparticles of ~ 700 nm at the small- 
est nanoparticle density of 2 /im -2 . Thus, we rule out 
the possibility of finite size effects in the following analy- 
ses. The samples were introduced into a vacuum chamber 
with a base pressure of ~ 10 -7 Torr and annealed at ~ 
500 °C for ~ 5 hours in order to remove any adhesive tape 
residue and to achieve equilibrium structures. After the 
annealing procedure, we observed surfaces of graphene 
flakes of various thicknesses in air using AFM in tapping 
mode with silicon cantilevers with a nominal tip radius 
of < 10 nm (Olympus, OMCL- AC160TS) . 



III. RESULTS AND DISCUSSION 

A. Wrinkling of monolayer graphene 

Figure [I] shows typical AFM images of monolayer 
graphene supported on nanoparticles for various densities 
p np . At p np = 11 /im -2 (Fig. [IJi), graphene adheres con- 
formally to the substrate, as noted previously [17, 28H30], 
with predominantly isolated protrusions at the nanopar- 
ticle locations. At p np = 22 /im -2 (Fig. [TJd), some 
nanoparticle-induced protrusions are linked by wrinkles 
(we use the term "wrinkle" in accordance with the lit- 
erature on graphene). Additional wrinkles with one free 




w 



FIG. 2. Schematics of (a) a wrinkle formed between two 
nanoparticles with diameters d and (b) the wrinkle profile 
along the transverse direction as represented by shaded area 
in (a). 



termination are also observed. The wrinkles between pro- 
trusions sag. If the protrusions have comparable heights 
as we assume below for simplicity, the wrinkle sags in 
the middle, while if the height difference is large, the 
wrinkle sags asymmetrically toward the protrusion of 
smaller height. With further increase in nanoparticle 
density, the wrinkles connecting the protrusions prolif- 
erate (Fig. pj) , and ultimately a wrinkle network spans 
the sample (Figs. [T]i and e). These observations indicate 
the presence of a critical distance X c between nanopar- 
ticles, below which wrinkling is induced. 

We now analyze the critical nanoparticle separation 
X c . The ridge running along the wrinkle between two 
nanoparticles with diameters d separated by X follows a 
catenary-like profile with a deflection ((x) as shown in 
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Fig. [2^i. Additionally, as represented in Fig. [2Jd, the pro- 
file of the ridge along the transverse (y) direction can be 
characterized with the dihedral angle and the curvature 
radius Co(x)~ 1 . The contour of the wrinkle results from 
the balance between elasticity and adhesion. A wrinkle of 
length X much larger than its average width W av g costs 
adhesion energy ~ TXW avg , where T is the graphene- 
Si02 adhesion energy per area. The deflection of elastic 
sheets between two protrusions creates a stretching strain 
6 ~ (Co /^avg) 2 along the crease of the ridge in a region of 
width ~ Co [31], where Co = C(0) is the maximum deflec- 
tion. Balancing the stretching energy ~ i^D^Co^ 2 with 
the adhesion energy, and assuming that Wavg ~ d, one 
finds that Co ~ X^ b d 1 ^(Y/E 2D ) 1 ^. Precluding wrin- 
kles with deflections larger than one finds a maximum 
wrinkle X c ~ d(E 2D /T) 1 ^. 

We next present a detailed elastic analysis of the wrin- 
kle shape including bending energy, which provides an 
expression for the full deflection profile ((x) and recov- 
ers the scaling law for X c . Assuming the opening angle 
6 is independent of x as validated in Ref. [3lJ we find 
the width of the deformed region w(x) = (tt — 0)Co(x)~ 1 
and the deflection CO) = [l/sin(0/2) - l]C o 0r) _1 within 
the effective one-dimensional model. Furthermore, we as- 
sume that the stretching strain in the y direction is irrel- 
evant according to Ref. [3TJ Then, the stretching strain 
is also given in one dimension (in the x direction) by 
e* = [1 + (dx() 2 } 1/2 ~ 1 ~ {d x Q 2 /2. We find the stretch- 
ing energy E s and the bending energy E^: 
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The adhesion energy cost is proportional to the area of 
the substrate uncovered by the membrane: 
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FIG. 3. (a) Profile of the wrinkle along the white dotted 
line in the AFM image shown in the inset. The scale bar 
in the inset is 50 nm. The solid red lines are theoretical 
expectations, (b) The maximum deflection Co as a function of 
the wrinkle length. The error bar indicates the uncertainty of 
Co due to the height difference between the protrusions. We 
choose wrinkles formed between protrusions where the height 
differences are less than 10 %. The area shaded in red is the 
theoretical prediction for scaling of Co with X for F — 0.6 — 2.8 
eV/nm 2 . The inset is a typical AFM image of the wrinkles 
formed between the two protrusions. The scale bar is 20 nm. 



E T = T J dxW, (3) 

where W is the base of the wrinkle profile as illustrated 
in Fig. [2}}. In addition, we take into account bending and 
adhesion at the foot of the wrinkle, which cost bending 
energy Ey and adhesion energy Ey> (see Appendix A for 
details), but these turn out to be negligible. 

At equilibrium, we expect SE tot /S( = 0, where E tot = 
E s + E b + E r + Ey + E v < , leading to a differential equa- 
tion for deflection ( s . By solving the differential equation 
under two boundary conditions £(±X/2) = 0, we find 
the deflection on both sides of the center of the wrinkle 
(see Appendix A): 

r (27k V /6 [ 1 l 1/3 /X \ 2/3 . 



In Fig. [3^i, we show the line profile along a wrinkle formed 
between two protrusions. The observed deflection is well 
fitted by the theoretical prediction (±(x) ~ (X/2=F x) 2 / 3 
with a prefactor of 0.32 nm 1 / 3 . Then, using d — 7.4 ± 
2.2 nm, E 2D = 2.12 x 10 3 eV/nm 2 [TO], k = 1 eV [9], 
and T = 0.6 - 2.8 eV/nm 2 [T71 [28l [32], we minimize 
the total energy E tot numerically with respect to for 
given X. Using Eq. [4j we find the maximum deflection 
Co = C(0) as a function of X as illustrated in Fig. [3J3- The 
maximum deflection Co monotonically increases with X, 
in good agreement with the observations. The theoretical 
model for a deflection is based on the assumption that a 
wrinkle is formed between two sharp peaks. The finite 
sizes of the protrusions may be a cause of decrease of 
the deflection below the theoretically expected range in 
Fig. [3)3. Furthermore, we attribute the most likely source 
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FIG. 4. The density of wrinkles p w and the mean num- 
ber of wrinkles per protrusion p w /pn P (inset) as functions 
of nanoparticle density p np . Each arrow corresponds to the 
AFM images shown in Figs. la-e. The solid red lines are fits 
described in text. 



of uncertainty to the observed dispersion in nanoparticle 
sizes. 

Since a wrinkle is geometrically suppressed if £(0) > d, 
the maximum length of the wrinkle is determined by a 
condition that £(0) = d. From Eq. [4j we find the maxi- 
mum length X c — 104 — 65 nm along with 6 = 35° — 14° 
for the adhesion energy T = 0.6 — 2.8 eV/nm 2 , respec- 
tively, in rough agreement with the observed maximum 
wrinkle length of ~ 200 nm [27] . Neglecting the contribu- 
tions of Ey and Ey> , X c slightly decreases to 96 — 62 nm. 
The discrepancy between the theoretical predictions and 
the observations is likely due to the fluctuations in the 
nanoparticle sizes which strongly influence the wrinkle 
length X c (X c increases in proportion to d). In Appendix 
B, we discuss the link between the detailed elastic model 
and the scaling analysis. 



B. Random wrinkling model 

Figure [4] shows the density of wrinkles p w as a function 
of the nanoparticle density p np and the number of wrin- 
kles per protrusion p w /p np as a function of p np (inset). 
We find p w is almost zero below p np ^25 pm~ 2 (arrow 
b) and then begins to increase rather linearly with p np 
above p np « 50 pm~ 2 (arrow c). We now analyze the be- 
havior of the wrinkle density p w versus the nanoparticle 
density p np within a simple model. 

We first consider nanoparticles placed at random on 
the substrate [27]. Then wrinkles are placed with a 
probability Q w between neighboring nanoparticles sep- 
arated by less than a cutoff length X c . The probability 
Q w is expected to encompass not only the true prob- 
ability to make a wrinkle between nanoparticles, but 
also all the information about wrinkle-orientation cor- 
relations, selecting only a fraction of all possible wrin- 



kles. Since nanoparticles with more than three con- 
nected wrinkles are scarcely observed [27 , we set three 
as the maximum number of wrinkles. Employing the 
probability density for a nanoparticle to have the z-th 
nearest nanoparticle (i — 1,2, and 3) at a distance r, 
Pi(r) = 2 (7rp np y +1 r 2z+1 exp (-7rp np r 2 ) ji\ [33, we find 
the density of wrinkles; 



2~~ 

Pnp^"u 



dr Pi(r) 



-nX 2 cPnp ( 2 + -irX 2 cPnp ) e"^"^ 



3 (l - e - '""*^) 



(5) 



The factor of 1/2 in p w compensates for the double- 
counting of each wrinkle (i.e. from the particles at each 
end). In the small nanoparticle density limit p np <C X~ 2 , 
the density of wrinkles is p w = (l/2)^ w 7rX 2 p 2 p , while in 
the large density limit p np ^> X~ 2 , each nanoparticle has 
at least three neighboring nanoparticles within distance 
X c , leading to p w = (3/2)Q w p np . The red solid lines 
in Fig. [4] are fits to Eq. [5] with & w = 0.54 and X c = 120 
nm. The cutoff length is consistent with the observations. 
Furthermore, the agreement with X c predicted from the 
above elastic analysis is good. The model indicates a sig- 
nificant increase of the wrinkle density for the nanopar- 
ticle density larger than (7rX c ) -2 , but also suggests that 
p w does not exhibit any singularity, i.e. wrinkling is a 
crossover phenomenon rather than a sharp transition. 



C. Percolation transition in the wrinkle network 

With increasing p npi the connectivity of the wrinkle 
network increases, and we find a percolation transition 
at a threshold density p c (of order X~ 2 ) at which the 
wrinkle network spans the system, as shown in the inset 
of Fig. [5^. The expansion of the network via wrinkling 
is a purely two-dimensional (2D) phenomenon. Thus, we 
analyze this behavior using 2D percolation theory [34] : In 
Fig. [5^i, we plot the probability P that a given protrusion 
belongs to the percolating cluster spanning a region of 
size LxL, where L ranges from 1 to 3 pm. Also plotted 
is the prediction from 2D percolation theory: P ~ (p np — 
p c Y for p np > p c with p c = 87.5 pm~ 2 as determined 
below and the "magnetization" exponent f3 = 5/36 [34] , 
which reproduces the observations well. In Fig. (5}d, we 
show the probability II that a cluster connects opposite 
sides of a region of size L x L (L = 0.5, 1, 2, and 3 /im). 
For an infinite system, II = 1 for p np > p c , while II = 
for p np < p c [34 . Indeed II displays a sharp transition 
around p np = 87.5 pm~ 2 for L = 3 /im, indicating p c is 
in that vicinity. 

Next, we probe the width A of the transition region, 
which is expected to scale as L -1 /^, where v — 4/3 is 
the correlation- length exponent [34]. We define A as the 
difference in density for II = 0.9 and II = 0.1 in Fig. [5b. 
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FIG. 5. (a) P as a function of nanoparticle density for L = 1, 2, and 3 /im. The inset is an AFM image (lxl /im 2 ) of graphene 
on the Si02 nanoparticles with a density of 57 /zm 2 , showing the percolating cluster highlighted by the dashed curve, (b) II as 
a function of the density of nanoparticles for L = 0.5, 1, 2, and 3 /im. Points for L = 0.5, 1 and 2 /mi represent averages in a 
bin of 10 /im -2 . The inset is a plot of log A as a function of logL; the red line indicates a best-fit power exponent of —1.0. (c) 
The mean finite cluster size S as a function of the density of nanoparticles (points represent averages in a bin of 2 /im -2 ). The 
red dashed line is the theoretical expectation (described in text). 



The inset of Fig. [SJd shows that the data are well-fitted 
with v — 1.0 ±0.3, consistent with the theoretical expec- 
tation. 

Finally, we plot in Fig. the mean size S of the clus- 
ters (excluding the percolation cluster) as a function of 
p np with the theoretical prediction £ ~ \p np — p c |~ 7 , 
where 7 = 43/18 is the "susceptibility" exponent [34] . 
Some Monte Carlo simulations predict a much larger 
prefactor for p np < p c (e.g., a critical amplitude ratio 
of 50 ± 26 for a continuum model [35]), in reasonable 
agreement with the observed ratio of ~ 30. 

Thus, all measurements strongly support the existence 
of a 2D percolation transition at a critical nanoparticle 
density p c « 87.5 /im -2 . Since the only length scale is 
X c , we obtain a universal number (i.e. independent of 
model parameters such as T, E2D, or d) characterizing 
the wrinkle percolation transition: p c X 2 « 0.9. In con- 
trast, simple continuum percolation of penetrable discs 
of diameter X c leads to p c X^ w 2.9 [35]. This difference 
is a consequence of the unique structure of the wrinkle 
network (e.g., not more than three wrinkles merging at 
a given nanoparticle [27]). 

D. Delamination of graphene multilayers 

Finally, we investigate morphological transitions 
which occur in multilayer graphene, using the same 
nanoparticle-templated substrates. The capability of 
multiple layers of graphene to mechanically screen an as- 
perity recalls the use of multiple mattresses in an attempt 
to hide the presence of a pea in the fairy tale "Princess 
and the Pea". Figures [6^1- f show typical AFM images of 
mono- and multi- layer graphene supported on nanopar- 
ticles of density 160 ± 24 /im -2 . The thicker graphene is 
partially suspended over the nanoparticles, as schemat- 



ically shown in the insets of Figs. [6]i and [6]f, with the 
delaminated area increasing with layer number n. 

The phase image in Fig. [7] demonstrates that the 
suspended graphene is indeed supported by isolated 
nanoparticles; the mechanical response of the graphene to 
the AFM tip allows the detection of the hidden nanopar- 
ticle "peas" under the flat graphene "mattress". The 
phase image records the varying phase angle of the (os- 
cillating) AFM cantilever as it interacts with an inhomo- 
geneous sample surface. The phase angle increases with 
increasing local sample stiffness [36]. Figure [7] shows that 
the phase image of 4-layer graphene discriminates be- 
tween rigid supported regions (larger phase) and flexible 
suspended regions (smaller phase). The high, flat regions 
in the topograph show small, roughly circular regions 
of large phase indicating the locations of the nanoparti- 
cles (arrows) which support the surrounding suspended 
graphene (small phase). 

Figure [6^ shows the areal fraction <j> of graphene in 
contact with the substrate and the characteristic length 
I of the delaminated regions, as functions of n. As n 
increases, a first transition occurs around n = 10, where I 
increases rapidly (Fig. [6]i; partial delamination); second, 
(j) decreases and becomes negligibly small above n ~ 15 
(Fig. [6]f; complete delamination). 

Surface-roughness-induced delamination of graphene 
has recently been studied theoretically [Til [El US] and 
experimentally [20l [2TJ [37]. Models assume the elastic 
energy is dominated either by bending [11] or stretch- 
ing [15]. Here we consider each regime, and assume that 
the adhesion energy between Si02 and n-layer graphene 
T n is independent of n for n > 1 and has the value 1.9 
eV/nm 2 [32]. In the bending-dominated model [11 , un- 
binding is controlled by a single dimensionless parameter 
a = (2r n //^ n ) 1 / 4 /[27r(p np (i) 1 / 2 ], where n n is the bending 
rigidity of n-layer graphene for n > 1. Without inter- 
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FIG. 6. Typical AFM images of (a) monolayer, (b) trilayer, (c) 7-layer, (d) 10-layer, (e) 14-layer, and (f) 18-layer graphene 
on Si02 with nanoparticle density 160 =b 24 /im" 2 . The scale bar in each image is 400 nm. The insets in (a), (d), and (f) are 
corresponding schematics of graphene films on nanoparticles to the AFM images, (g) The fractional area cj) where graphene 
conforms to the substrate and the characteristic length I of the delaminated domains, as functions of the number n of graphene 
layers. The dashed curves are the theoretical expectations for cj) (orange) and / (blue) (described in text). 
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FIG. 7. Typical AFM (a) height and (b) phase images (lxl 
/im 2 ) of 4- layer graphene supported on the nanoparticles. The 
scale bar is 200 nm. (c) Line profiles of the height and phase 
along the dashed red and blue lines shown in (a) and (b), 
respectively. The arrows correspond to those in (b), showing 
the locations of nanoparticles beneath graphene. 



layer sliding [38 , continuum plate elasticity [39 gives 
n n = Et 3 n 3 /[12(1 - v 2 g % where t = 0.335 nm is the in- 
terlayer spacing, E w 0.96 TPa is the Young's modulus, 
and v g ~ 0.165 is Poisson's ratio of monolayer graphene 
[TO] . The threshold for partial unbinding is predicted at 



0.8 < a < 1.3, or 3 < n < 6, with complete unbinding 
at 0.55 < a < 0.75, corresponding to 7 < n < 10 [TT] . 
Thus, the bending-dominated model underestimates the 
critical value of n for unbinding, indicating that it overes- 
timates the bending elastic energy. The one-dimensional 
character of the bending model limits its ability to make 
quantitative predictions. Furthermore, given the small 
radii of curvature in our experiment, the bending energy 
might well be reduced by partial interlayer sliding. Per- 
fect sliding would give n n = uk, leading to an unbinding 
threshold for n a hundredfold larger; hence, interlayer 
sliding is extremely effective in relieving bending stress. 

We then develop a stretching-dominated model, us- 
ing Schwerin's solution for a membrane pushed by a 
point force [40]. The solution gives the diameter of the 
detachment zone in n-layer graphene around a protru- 
sion: 2R « d(4nE 2D /3T n ) 1 / 4 (Refs. [15] and HQ see also 
Appendix C), where nE2D is the tensile rigidity of n- 
layer graphene. Then, we notice that the detached area 
around each protrusion is ttR 2 and, furthermore, we as- 
sume that the detached areas produced by the wrinkles 
are negligible. Therefore, the typical length of the de- 
laminated regions / is estimated to be 2R. The adhered 
area fraction <j) is equivalent to the probability to have 
no nanoparticle in a domain of an area of ttR 2 , leading 
to (j) = exp(— i:R 2 p np ) (see Appendix D for details). As 
shown in Fig. [6^, these predictions reproduce well the 
observed thickness dependence of <j) and I below n ~ 10, 
indicating that the stretching-dominated model for iso- 
lated protrusions accurately describes the small-n limit 
where p np <C l~ 2 . However, / increases and (j) decreases 
much more rapidly than these predictions for n > 10, 
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indicating collective effects become important. In order 
to understand the collective delamination in the high- 
nanoparticle-density regime p np > l~ 2 , we may need to 
solve full elastic membrane equations, i.e., the Foppl-von 
Karman equations [39] allowing for multiple nanoparti- 
cles. 



IV. PSEUDOMAGNETIC FIELDS IN 
WRINKLED GRAPHENE 

Finally, while a complete treatment is beyond the 
scope of this paper, we briefly discuss the potential im- 
pact of wrinkles on the electronic properties of graphene. 
Strain in graphene produces a vector potential in the 
electronic Hamiltonian, with strain gradients resulting in 
effective magnetic fields [23] . Such pseudomagnetic fields 
generated by strain gradients have been proposed as the 
basis of engineering graphene's electronic properties [241 - 
[26]. 

We first evaluate the pseudomagnetic field generated 
by strain in the absence of wrinkling, corresponding to 
the case of small thickness or small nanoparticle density. 
In this case the elastic behavior of graphene on nanoparti- 
cles is predominantly determined by stretching, resulting 
in significant strain. The radial strain e r and the circum- 
ferential strain both scale as ~ (Td 2 / E2d) 1 ^ r -2 / 3 
at < r < R, although the radial strain is more than 
5 times as large as the circumferential strain (see Ap- 
pendix C and Ref. [42]) . The gauge fields induced by 
the axisymmetric strain can be written as (A^A^) « 
($o/3/a)(e r - e^)(cos3(/?,sin3(/?), where $ = 10~ 15 Wb 
is the flux quantum, /3 ~ 2 is the change in the hopping 
amplitude between the neighboring atomic sites due to 
the lattice deformation [43], a = 0.142 nm is the lattice 
constant, and cp is the azimuthal angle with cp = in the 
zig-zag direction [25 . Then, the strain-induced pseudo- 
magnetic field is given by £? e ff = d<pA r /r — d r A^ — A^/r. 

Thus, we find £ eff - ($ /V a ) {^d 2 /E 2D ) 1/S r~ 5/3 sin3<p. 
The strain in graphene on an isolated nanoparticle in- 
duces three-fold symmetric pseudomagnetic field profiles 
with maximum fields along the arm-chair directions. The 
pseudomagnetic field pattern is similar to a recent exper- 
iment [44], in which suspended graphene was deformed 
by a sharp tip and pseudomagnetic fields were found to 
confine electrons to quantum dots with charging energies 
and level spacings both of order 10 meV. The divergence 
of the strain at r = is cut off by the finite radius of 
the nanoparticles; thus, we may expect that maximum 
pseudomagnetic field appears at a radius comparable to 
the nanoparticle radius. We therefore estimate that the 
maximum pseudomagnetic field £? e ff is of order 300 T for 
d = 7.4 nm and r = d/2, significantly greater than in 
Ref. [44j suggesting that the impact on electronic proper- 
ties may be even more profound. 

Next we evaluate strain and strain-induced pseudo- 
magnetic fields in a wrinkle. Strain along a wrinkle is 
given by e x ~ (d x () 2 . Then, using Eq. [ZJ we find 



the strain distribution e x ~ / 'E 2 dY " /a \X / '2 =p x)" 2/3 . 
The pseudomagnetic field is estimated to be B e & « 
$of3t x /(aW) [15] [23], where W is the typical wrinkle 
width. In the strong adhesion limit d^T/n) 1 / 2 ^> 1, the 
wrinkle width W can be estimated to be ~ (hz/2T) 1 ^ 2 « 1 
nm (see Appendix B). Thus, we find in the middle of a 
wrinkle the pseudomagnetic field has a broad minimum 
on the order of 10 T for X = 100 nm. 10 T is a large 
magnetic field compared to the disorder strength 1 //i ~ 1 
T in typical graphene samples (/i being the electron mo- 
bility) and corresponds to an energy difference between 
0th and 1st Landau levels of ~ 1300 K. Hence, we expect 
pseudomagnetic field effects due to wrinkles in graphene 
to be significant. 

The pseudomagnetic field near particles in the wrin- 
kled case will generally be more complicated, depending 
on the number of wrinkles terminating on the particle 
and their direction with respect to each other and the 
lattice. However, qualitatively we expect that since wrin- 
kling reduces the in-plane strain around the nanoparti- 
cles, the pseudomagnetic field is also reduced. Recent 
molecular dynamics simulation results [45] relevant to 
experiments by Tomori et al. [46 and our experiments 
of graphene on nanoparticles have indeed revealed that 
when nano-scale pillars supporting graphene are located 
far away from each other, graphene is detached only 
around the pillars and three-fold symmetric pseudomag- 
netic fields are induced around each pillar, while with 
decreasing distance between the pillars, graphene delam- 
inates in regions between the pillars, resulting in com- 
plicated pseudomagnetic field profiles. Our observations 
of wrinkling and delamination combined with theoreti- 
cal analysis based on a continuum elastic model can be 
used to place limits on strain distributions and, thus, 
on pseudomagnetic field maxima realizable in monolayer 
graphene through adhesion to patterned surfaces. 



V. CONCLUSIONS 

In conclusion, combining experiments with monolayer 
and multilayer graphene, we have obtained a global 
picture of the structural evolution of graphene mem- 
branes on surfaces of varying roughness. With increas- 
ing nanoparticle density (or graphene thickness), this 
evolution proceeds in five stages: conformal adhesion, 
wrinkling, wrinkle percolation, partial delamination, and 
complete delamination. The results can be used to place 
upper limits on the magnitude of pseudomagnetic fields 
generated in graphene by adhesion to patterned surfaces. 
Finally, the wrinkling and delamination are not specific 
to graphene; they should be a general feature of soft 
membranes adhered to rough surfaces, with implications 
for systems ranging from cell walls to fabrics. 
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Appendix A: Detailed geometric model for the 
wrinkle shape 



with the two boundary conditions ((±X/2) — 0. We an- 
ticipate £ to be symmetric with respect to x, so that d x ( 
should vanish at x = 0. However, Eq. A6 indicates that 
if d x ( vanishes at x = 0, either £ or d zz ( should diverge. 
Since the solution with diverging ( is physically incon- 
ceivable, we conclude that d zz ( should diverge. This in- 
dicates a discontinuity of the slope at x = 0. Physically 
this singularity would be regularized at small scales ei- 
ther by bending along the x direction or by stretching 
along the y direction. These contributions are expected 
to be small. As a result, we obtain a simple solution on 
both sides of the center of the wrinkle as described by 
Eq.|U 



The stretching energy E s , the bending energy E^, and 
the adhesion energy Er given by Eqs. [l][3]are reduced to 



E s = 2 ^ J dxw(x)e 2 : 

E2D / „x r 1 



(7T-0) 



sin(0/2) 



dx((d x () 4 



(Al) 



1 



sin(0/2) 



dxC 1 (A2) 



E r = T J dxW = 2rXdtan(<9/2). (A3) 

The combination between the bending and the adhe- 
sion at the foot of the wrinkle forces an equilibrium cur- 
vature at the contact line, C eq = (2r/^) 1/2 [47]. Then, 
assuming a simple geometry where a constant curvature 
region matches the constant slope region as shown in 
Fig. [2}d, we find the bending energy Ey and the adhe- 
sion energy Ey> of the curved zone: 



Ey = 2- I dx I dyC 2 eq =X 



1/2 



(tt - 0) (A4) 



Er> = 2TXC~ 1 tan/3 = X(2Yk) 1/2 tan 

\ 1 / 

(A5) 

where 2/3 is the angle of the curved region as shown in 
Fig. g>. 

Minimization of the total energy E tot = E s + E\> + 
E r + Ey + E T > with respect to £ leads to the following 
differential equation: 



C 2 [3(^C) 4 +i2C(^C) 2 ^C] 



4k 



E' 



2D 



1 

sin(0/2) 



= 
(A6) 



Appendix B: Scaling analysis for X c 

Here we show that X c scales as X c ~ d(E2D /Y) 1 ^ , 
analogous to scaling for the diameter detachment zones 
surrounding a local protuberance [15] . We design a scal- 
ing analysis, neglecting Ey> and Ey . The total energy is 
of the form 

E tot = ^E\>«X^[h{e) + TdK-^E-^X^h{9)}, 

(Bl) 

where 



sin(6>/2) 



1 



2/3 



and 



/ 2 (0) = 2tan-. 



(B2) 



(B3) 



The minimization of Eq. |B1| with respect to leads to 

deh{e) + TdK- 5 / 6 E-^ 6 X 2 / 3 def 2 (0) = (B4) 

and as a consequence 

e = h(Yd K -^E-^X 2 l*). (B5) 

The critical spacing X c is given by a condition that 
£(0) = d. Hence, using Eq. [4| we find 



o3/2 / p \ 1/4 

33/4 



d 3/2 ( E2D 



E2D 



1/4 
1/4 



1 



sin(0(X c )/2) 



-1/2 



h(fs(Td^E^ 6 X^)l 



where 



U(0) 



1 

_sin(0/2) 



-1/2 



(B6) 



(B7) 



One can check by substitution that 

1/4 



Xc = d3/2 (^T 



h(Td 2 /K) 



(B8) 



with f 6 (u) = h(uh(u) 2 t*) and f 5 (u) = U(f 3 (u)). We 
now define the elastic thickness h e \ = {^/E 2 d) 1 ' 2 an d the 
equilibrium contact curvature C eq = (2r/ ft) 1 / 2 . Letting 
/V^) = /q(u 2 /2)^ we rewrite JT C as 



d 3/2 

-jj^f7(C eq d), 
hi 



(B9) 




FIG. 8. Schematic of graphene supported on a single nanopar- 
ticle. The detachment length is R. 



which is the general scaling form of the solution. 

We now consider two asymptotic limits; the strong 
adhesion limit C eq d ^> 1 and the weak adhesion limit 
C eq d <C 1. In the strong adhesion limit, the opening an- 
gle of the wrinkle goes to zero. Then, one has fi(6) ~ 
6>- 2 / 3 and f 2 (0) ~ 0. Hence f 3 (u) ~ u~ 3 ^. Since / 4 (0) - 
Pl\ one has f 5 (u) = h{h{u)) ~ [/ 3 H] 1/2 ~ ^ 3/1 °, 



and finally fe(u) 



_1 / 4 . Therefore one has: 



X, 



^3/2 



(C eq d) 



-1/2 



and 



(c eg d)- 



d(E 2D /r) 1 / 4 (BIO) 



(Bll) 



which also confirms that the small-(9 limit corresponds to 
the large-Cegd limit (strong adhesion limit). 

Alternatively, setting = ir — r9, we redo the above 
scaling analysis in the weak adhesion limit C eq d <C 1 
with the argument of /i, f 2l and f± being 6 instead of 0. 
Then A(6) - 9 7 / 3 , and / 2 (6) - 9" 1 , so that f 3 (u) - 
ii 3 / 10 . Also / 4 (6) - 6- 1 , leading to f 5 (u) ~ u~ 3 ^ and 
fe(u) ~ u~ x l^. Consequently, we obtain once again 



d 3 ' 2 
iM 2 



(C eq d) 



-1/2 



d(E 2D /T) 



1/4 



and 







(C eq d) 1 ^. 



(B12) 



(B13) 



This solution is consistent with the weak adhesion limit 
because 6 < 1 implies C eq d <C 1. 



Appendix C: Stretching of graphene on a 
nanoparticle 

The typical extent of the detachment zone of graphene 
caused by the presence of a nanoparticle is discussed here. 
In the simplest model, we use the power-law solution of 
Schwerin [40], in the regime dominated by stretching. 
In this regime, the bending rigidity only contributes as 



boundary layer effects at the attachment lines [42] . Con- 
sequently, the Schwerin solution does not match tangen- 
tially to the substrate and the nanoparticle (the match- 
ing being forced at a smaller scale controlled by bend- 
ing rigidity). Assuming that the nanoparticle diame- 
ter d is much smaller than the radius R of the detach- 
ment zone (see Fig. [8J, we obtain the angle of rota- 
tion rj(r) = [8F/ (9tv E 2 Dr)] 1 ^ 3 and the vertical distance 
from the substrate Z = [3R 2 F/(7rE 2D )}^ 3 (see Fig. § 
and Ref. |42|) , where F is the force on the apex. How- 
ever, the above solution does not match the boundary 
conditions, and a better approximate numerical solution 
is: Z = g(v g )(R 2 F/E 2D ) 1 / 3 , where g(v g ) = 1.0491 - 
0.1462^-0.15827^ [42 . For graphene, v g =0.165 [TO], 
and g{v g ) « 1.029 is very close to (3/tt) 1 / 3 « 0.984 so 
that we can safely use the direct Schwerin solution. 

The elastic stretching energy can be calculated from 
a gedanken experiment, where the height Z is increased 
with constant R: 



E(Z) 



f 

Jo 



7rE< 



dzF(z) 



2D 



3R 2 jo 
irE 2D Z 4 
12R 2 



dzz 6 



(CI) 



Assuming that the apex height is equal to the diameter d 
of the nanoparticle, one has Z = d, and the total energy 
reads 



E(d) 



nE 2D d A 
12R 2 



TirR 



(C2) 



Minimizing the total energy with respect to i?, one finds 
2R/d = (4E2D/3T) 1 / 4 as suggested in Ref. ED 

Assuming that a pre-existing strain is negligibly small, 
one can find radial and circumference components of 
strain present in an isolated graphene protrusion at < 
r < R: 



e r (r) 



— ( G 

E 2D \ r 



dG 

— COS 77 + H sin 77 — v a — — 
dr 



z-u g f 4r \ 1/3 fd\ 

4 \9E 2D J \rj 



2/3 



(C3) 
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E2D 
1 - 3i/ £ 



dG 
dr 



G 



COS T] 

1/3 



9 



v g H sin r] 

2/3 



(C4) 



4 V 9 ^2D/ 

where G = F/{2nr]) = 9^2 j d^ 2 ^/16 is the stress variable 
in the radial direction, H = F/(27rr) = 9E 2 dV 3 /^ is the 
vertical sheet stress, and we assume 7?~0at0<r<i? 

m. 



Q(A)-dA p np Q(A), leading to dQ(A)/dA = -p np Q(A). 
Using the normalization constraint 



lim Q(A) 



1, 



(Dl) 



Appendix D: Area probability distribution 

The probability Q(A) that a nanoparticle has no 
nanoparticles in a domain of area A obeys Q(A + dA) = we find Q(A) = exp(— p np A). 
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